####Model results####
##a) demographic plurality (table X35, models A123-A126)
r4b_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_largest_g == 1))
r4b_m1g_t1_sc_cse <- data.frame(cluster.se(r4b_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_largest_g == 1)$cowcode))[,2])
r4b_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_largest_g == 0))
r4b_m1g_t1_sf_cse <- data.frame(cluster.se(r4b_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_largest_g == 0)$cowcode))[,2])
r4b_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_largest == 1))
r4b_m1d_t1_sb_cse <- data.frame(cluster.se(r4b_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_largest == 1)$cowcode))[, 2])
r4b_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_largest == 1))
r4b_m2d_t1_sb_cse <- data.frame(cluster.se(r4b_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_largest == 1)$cowcode))[, 2])
stargazer(r4b_m1g_t1_sc, r4b_m1g_t1_sf, r4b_m1d_t1_sb, r4b_m2d_t1_sb, se=c(r4b_m1g_t1_sc_cse, r4b_m1g_t1_sf_cse, r4b_m1d_t1_sb_cse, r4b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A123)", "Min. (model A124)","Maj./min. dyad (model A125)","Maj./min. dyad (model A126)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X35. Additional results: alternative definition of second-order majority (demographic plurality).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex35.txt")
stargazer(r4b_m1g_t1_sc, r4b_m1g_t1_sf, r4b_m1d_t1_sb, r4b_m2d_t1_sb, se=c(r4b_m1g_t1_sc_cse, r4b_m1g_t1_sf_cse, r4b_m1d_t1_sb_cse, r4b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A123)", "Min. (model A124)","Maj./min. dyad (model A125)","Maj./min. dyad (model A126)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X35. Additional results: alternative definition of second-order majority (demographic plurality).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex35.html")
##b) demographic majority (table X36, models A127-A130)
r4_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_majority_g == 1))
r4_m1g_t1_sc_cse <- data.frame(cluster.se(r4_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_majority_g == 1)$cowcode))[,2])
r4_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_majority_g == 0))
r4_m1g_t1_sf_cse <- data.frame(cluster.se(r4_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_majority_g == 0)$cowcode))[,2])
r4_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_majority == 1))
r4_m1d_t1_sb_cse <- data.frame(cluster.se(r4_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_majority == 1)$cowcode))[, 2])
r4_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_majority == 1))
r4_m2d_t1_sb_cse <- data.frame(cluster.se(r4_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_majority == 1)$cowcode))[, 2])
stargazer(r4_m1g_t1_sc, r4_m1g_t1_sf, r4_m1d_t1_sb, r4_m2d_t1_sb, se=c(r4_m1g_t1_sc_cse, r4_m1g_t1_sf_cse, r4_m1d_t1_sb_cse, r4_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A127)", "Min. (model A128)","Maj./min. dyad (model A129)","Maj./min. dyad (model A130)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X36. Additional results: alternative definition of second-order majority (demographic majority).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex36.txt")
stargazer(r4_m1g_t1_sc, r4_m1g_t1_sf, r4_m1d_t1_sb, r4_m2d_t1_sb, se=c(r4_m1g_t1_sc_cse, r4_m1g_t1_sf_cse, r4_m1d_t1_sb_cse, r4_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A127)", "Min. (model A128)","Maj./min. dyad (model A129)","Maj./min. dyad (model A130)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X36. Additional results: alternative definition of second-order majority (demographic majority).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex36.html")

####Figure A19####
##a) demographic plurality (civil)
me_r4b_m1g_t1_sc <- margins_summary(r4b_m1g_t1_sc, variables = "sa_territory_t", vcov = cluster.vcov(r4b_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_largest_g == 1)$cowcode)))
me_r4b_m1g_t1_sc$term <- "second-order maj."
me_r4b_m1g_t1_sf <- margins_summary(r4b_m1g_t1_sf, variables = "sa_territory_t", vcov = cluster.vcov(r4b_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_largest_g == 0)$cowcode)))
me_r4b_m1g_t1_sf$term <- "second-order min."
me_r4b_m1g <- rbind.fill(me_r4b_m1g_t1_sc, me_r4b_m1g_t1_sf)
me_r4b_m1g$term <- as.factor(me_r4b_m1g$term)
me_r4b_m1g$model <- "demographic plurality"
##b) demographic plurality (communal)
me_r4b_m1d_t1_sb <- margins_summary(r4b_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(r4b_m1d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & int_grp_rel_largest == 1)$cowcode)))
me_r4b_m1d_t1_sb$term <- "all"
me_r4b_m1d_t1_sb$model <- "demographic plurality"
me_r4b_m2d_t1_sb <- margins_summary(r4b_m2d_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r4b_m2d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & int_grp_rel_largest == 1)$cowcode)))
me_r4b_m2d_t1_sb$term <- ifelse(me_r4b_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r4b_m2d_t1_sb$model <- "demographic plurality"
##c) demographic majority (civil)
me_r4_m1g_t1_sc <- margins_summary(r4_m1g_t1_sc, variables = "sa_territory_t", vcov = cluster.vcov(r4_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_majority_g == 1)$cowcode)))
me_r4_m1g_t1_sc$term <- "second-order maj."
me_r4_m1g_t1_sf <- margins_summary(r4_m1g_t1_sf, variables = "sa_territory_t", vcov = cluster.vcov(r4_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_majority_g == 0)$cowcode)))
me_r4_m1g_t1_sf$term <- "second-order min."
me_r4_m1g <- rbind.fill(me_r4_m1g_t1_sc, me_r4_m1g_t1_sf)
me_r4_m1g$term <- as.factor(me_r4_m1g$term)
me_r4_m1g$model <- "demographic majority"
##d) demographic majority (communal)
me_r4_m1d_t1_sb <- margins_summary(r4_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(r4_m1d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & int_grp_rel_majority == 1)$cowcode)))
me_r4_m1d_t1_sb$term <- "all"
me_r4_m1d_t1_sb$model <- "demographic majority"
me_r4_m2d_t1_sb <- margins_summary(r4_m2d_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r4_m2d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & int_grp_rel_majority == 1)$cowcode)))
me_r4_m2d_t1_sb$term <- ifelse(me_r4_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r4_m2d_t1_sb$model <- "demographic majority"
##e) combined (civil)
me_demographic <- rbind.fill(me_r4b_m1g,me_r4_m1g)
me_demographic$estimate <- me_demographic$AME
me_demographic$conf.low <- me_demographic$lower
me_demographic$conf.high <- me_demographic$upper
me_demographic_plot <- dwplot(me_demographic,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("a) civil violence") +
  scale_color_manual(name = "model", values = c("#1b9e77","#d95f02","#7570b3"))+ 
  scale_linetype_manual(name = "model", values = c(3,2,1))+ 
  scale_shape_manual(name = "model", values = c(17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of civil violence") + ylab("group type")+ #+
  guides(shape = guide_legend(nrow=1,"model",reverse=TRUE), colour = guide_legend(nrow=1,"model",reverse=TRUE),linetype = guide_legend(nrow=1,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format(), breaks =  c(-0.015, -0.0075, 0))
##f) combined (communal)
me_demographic_d <- rbind.fill(me_r4b_m1d_t1_sb,me_r4b_m2d_t1_sb,me_r4_m1d_t1_sb,me_r4_m2d_t1_sb)
me_demographic_d <- subset(me_demographic_d, term != "other")
me_demographic_d$term <- as.factor(me_demographic_d$term)
me_demographic_d$estimate <- me_demographic_d$AME
me_demographic_d$conf.low <- me_demographic_d$lower
me_demographic_d$conf.high <- me_demographic_d$upper
me_demographic_d_plot <- dwplot(me_demographic_d,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("b) communal violence maj./min. dyad") +
  scale_color_manual(name = "model", values = c("#1b9e77","#d95f02","#7570b3"))+ 
  scale_linetype_manual(name = "model", values = c(3,2,1))+ 
  scale_shape_manual(name = "model", values = c(17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of communal violence") + ylab("dyad type") +
  scale_x_continuous(labels=scales::percent_format(), breaks = c(0, 0.0025, 0.005)) +
  guides(shape = guide_legend(nrow=1,"model",reverse=TRUE), colour = guide_legend(nrow=1,"model",reverse=TRUE),linetype = guide_legend(nrow=1,"model",reverse=TRUE))
##g) combined
figurea19 <- grid_arrange_shared_legend(me_demographic_plot, me_demographic_d_plot, ncol=2, nrow=1)
ggsave(figurea19, file='../figures/figurea19.pdf', width = 14, height = 4.5, units="cm",dpi=1000)